The present invention relates to the field of computer aided analysis of medical images. In particular, the present invention relates to a fast method for detecting lines in medical images.
Line detection is an important first step in many medical image processing algorithms. For example, line detection is an important early step of the algorithm disclosed in U.S. patent application Ser. No. 08/676,660, entitled xe2x80x9cMethod and Apparatus for Fast Detection of Spiculated Lesions in Digital Mammograms,xe2x80x9d filed Jul. 19, 1996, the contents of which are hereby incorporated by reference into the present application. Generally speaking, if the execution time of the line detection step can be shortened, then the execution time of the overall medical image processing algorithm employing that line detection step can be shortened.
In order to clearly illustrate the features and advantages of the preferred embodiments, the present disclosure will describe the line detection algorithms of both the prior art and the preferred embodiments in the context of the computer-assisted diagnosis system of U.S. patent application Ser. No. 08/676,660, supra. Importantly, however, the scope of the preferred embodiments is not so limited, the features and advantages of the preferred embodiments being applicable to a variety of image processing applications.
FIG. 1 shows steps performed by a computer-assisted diagnosis unit similar to that described in U.S. patent application Ser. No. 08/676,660, which is adapted to detect abnormal spiculations or lesions in digital mammograms. At step 102, an x-ray mammogram is scanned in and digitized into a digital mammogram. The digital mammogram may be, for example, a 4000xc3x975000 array of 12-bit gray scale pixel values. Such a digital mammogram would generally correspond to a typical 8xe2x80x3xc3x9710xe2x80x3 x-ray mammogram which has been digitized at 50 microns (0.05 mm) per pixel.
At step 104, which is generally an optional step, the digital mammogram image is locally averaged, using steps known in the art, down to a smaller size corresponding, for example, to a 200 micron (0.2 mm) spatial resolution. The resulting digital mammogram image that is processed by subsequent steps is thus approximately 1000xc3x971250 pixels. As is known in the art, a digital mammogram may be processed at different resolutions depending on the type of features being detected. If, for example, the scale of interest is near the order of magnitude 1 mm-10 mm, i.e., if lines on the order of 1 mm-10 mm are being detected, it is neither efficient nor necessary to process a full 50-micron (0.05 mm) resolution digital mammogram. Instead, the digital mammogram is processed at a lesser resolution such as 200 microns (0.2 mm) per pixel.
Generally speaking, it is to be appreciated that the advantages and features of the preferred embodiments disclosed infra are applicable independent of the size and spatial resolution of the digital mammogram image that is processed. Nevertheless, for clarity of disclosure, and without limiting the scope of the preferred embodiments, the digital mammogram images in the present disclosure, which will be denoted by the symbol I, will be Mxc3x97N arrays of 12-bit gray scale pixel values, with M and N having exemplary values of 1000 and 1250, respectively.
At step 106, line and direction detection is performed on the digital mammogram image I. At this step, an Mxc3x97N line image L(i,j) and an Mxc3x97N direction image xcex8max(i,j) are generated from the digital mammogram image I. The Mxc3x97N line image L(i,j) generated at step 106 comprises, for each pixel (i,j), line information in the form of a xe2x80x9c1xe2x80x9d if that pixel has a line passing through it, and a xe2x80x9c0xe2x80x9d otherwise. The Mxc3x97N direction image xcex8max(i,j) comprises, for those pixels (i,j) having a line image value of xe2x80x9c1xe2x80x9d, the estimated direction of the tangent to the line passing through the pixel (i,j). Alternatively, of course, the direction image xcex8max(i,j) may be adjusted by 90 degrees to correspond to the direction orthogonal to the line passing through the pixel (i,j).
At step 108, information in the line and direction images is processed for determining the locations and relative priority of spiculations in the digital mammogram image I. The early detection of spiculated lesions (xe2x80x9cspiculationsxe2x80x9d) in mammograms is of particular importance because a spiculated breast tumor has a relatively high probability of being malignant.
Finally, at step 110, the locations and relative priorities of suspicious spiculated lesions are output to a display device for viewing by a radiologist, thus drawing his or her attention to those areas. The radiologist may then closely examine the corresponding locations on the actual film x-ray mammogram. In this manner, the possibility of missed diagnosis due to human error is reduced.
One of the desired characteristics of a spiculation-detecting CAD system is high speed to allow processing of more x-ray mammograms in less time. As indicated by the steps of FIG. 1, if the execution time of the line and direction detection step 106 can be shortened, then the execution time of the overall mammogram spiculation detection algorithm can be shortened.
A first prior art method for generating line and direction images is generally disclosed in Gonzales and Wintz, Digital Image Processing (1987) at 333-34. This approach uses banks of filters, each filter being xe2x80x9ctunedxe2x80x9d to detect lines in a certain direction. Generally speaking, this xe2x80x9ctuningxe2x80x9d is achieved by making each filter kernel resemble a second-order directional derivative operator in that direction. Each filter kernel is separately convolved with the digital mammogram image I. Then, at each pixel (i,j), line orientation can be estimated by selecting the filter having the highest output at (i,j), and line magnitude may be estimated from that output and other filter outputs. The method can be generalized to lines having pixel widths greater than 1 in a multiscale representation shown in Daugman, xe2x80x9cComplete Discrete 2-D Gabor Transforms by Neural Networks for Image Analysis and Compression,xe2x80x9d IEEE Trans. ASSP, Vol. 36, pp. 1169-79 (1988).
The above filter-bank algorithms are computationally intensive, generally requiring a separate convolution operation for each orientation-selective filter in the filter bank. Additionally, the accuracy of the angle estimate depends on the number of filters in the filter bank, and thus there is an implicit tradeoff between the size of the filter bank (and thus total computational cost) and the accuracy of angle estimation.
A second prior art method of generating line and direction images is described in Karssemeijer, xe2x80x9cRecognition of Stellate Lesions in Digital Mammograms,xe2x80x9d Digital Mammography: Proceedings of the 2nd International Workshop on Digital Mammography, York, England, (Jul. 10-12, 1994) at 211-19, and in Karssemeijer, xe2x80x9cDetection of Stellate Distortions in Mammograms using Scale Space Operators,xe2x80x9d Information Processing in Medical Imaging 335-46 (Bizais et al., eds. 1995) at 335-46. A mathematical foundation for the Karssemeijer approach is found in Koenderink and van Doom, xe2x80x9cGeneric Neighborhood Operators,xe2x80x9d IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 14, No. 6 (June 1992) at 597-605. The contents of each of the above two Karssemeijer references and the above Koenderink reference are hereby incorporated by reference into the present application.
The Karssemeijer algorithm uses scale space theory to provide an accurate and more efficient method of line detection relative to the filter-bank method. More precisely, at a given level of spatial scale "sgr", Karssemeijer requires the convolution of only three kernels with the digital mammogram image I, the angle estimation at a pixel (i,j) then being derived as a trigonometric function of the three convolution results at (i,j) .
FIG. 2 shows steps for computing line and direction images in accordance with the Karssemeijer algorithm. At step 202, a spatial scale parameter "sgr" and a filter kernel size Nk are selected. The spatial scale parameter "sgr" dictates the width, in pixels, of a Gaussian kernel G(r,"sgr"), the equation for which is shown in Eq. (1):
G(r,"sgr")=(xc2xdxcfx80"sgr"2)exp(xe2x88x92r2/2"sgr"2)xe2x80x83xe2x80x83(1)
At step 202, the filter kernel size Nk, in pixels, is generally chosen to be large enough to contain the Gaussian kernel G(r,"sgr") in digital matrix form, it being understood that the function G(r,"sgr") becomes quite small very quickly. Generally speaking, the spatial scale parameter "sgr" corresponds, in an order-of-magnitude sense, to the size of the lines being detected. By way of example only, and not by way of limitation, for detecting 1 mm-10 mm lines in fibrous breast tissue in a 1000xc3x971250 digital mammogram at 200 micron (0.2 mm) resolution, the value of "sgr" may be selected as 1.5 pixels and the filter kernel size Nk may be selected as 11 pixels. For detecting different size lines or for greater certainty of results, the algorithm or portions thereof may be repeated using different values for "sgr" and the kernel size.
At step 204, three filter kernels K"sgr"(0), K"sgr"(60), and K"sgr"(120) are formed as the second order directional derivatives of the Gaussian kernel G(r,"sgr") at 0 degrees, 60 degrees, and 120 degrees, respectively. The three filter kernels K"sgr"(0), K"sgr"(60), and K"sgr"(120) are each of size Nk, each filter kernel thus containing Nkxc3x97Nk elements.
At step 206, the digital mammogram image I is separately convolved with each of the three filter kernels K"sgr"(0), K"sgr"(60), and K"sgr"(120) to produce three line operator functions W"sgr"(0), W"sgr"(60), and W"sgr"(120), respectively, as shown in Eq. (2):
W"sgr"(0)=I*K"sgr"(0)
W"sgr"(60)=I*K"sgr"(60)
W"sgr"(120)=I*K"sgr"(120)xe2x80x83xe2x80x83(2)
Each of the line operator functions W"sgr"(0), W"sgr"(60), and W"sgr"(120) is, of course, a two-dimensional array that is slightly larger than the original Mxc3x97N digital mammogram image array I due to the size Nk of the filter kernels.
Subsequent steps of the Karssemeijer algorithm are based on a relation shown in Koenderink, supra, which shows that an estimation function W"sgr"(xcex8) may be formed as a combination of the line operator functions W"sgr"(0), W"sgr"(60), and W"sgr"(120) as defined in equation (3):
W"sgr"(xcex8)=(⅓)(1+2 cos(2xcex8))W"sgr"(0)+(⅓)(1xe2x88x92cos(2xcex8)+(3)sin(2xcex8))W"sgr"(60)+(⅓)(1xe2x88x92cos(2xcex8)xe2x88x92(3)sin(2xcex8))W"sgr"(120)xe2x80x83xe2x80x83(3)
As indicated by the above definition, the estimation function W"sgr"(xcex8) is a function of three variables, the first two variables being pixel coordinates (i,j) and the third variable being an angle xcex8. For each pixel location (i,j), the estimation function W"sgr"(xcex8) represents a measurement of line strength at pixel (i,j) in the direction perpendicular to xcex8. According to the Karssemeijer method, an analytical expression for the extrema of W"sgr"(xcex8) with respect to xcex8, denoted xcex8min,max at a given pixel (i,j) is given by Eq. (4):
xcex8min,max=xc2xd[arc tan{(3)(W"sgr"(60)xe2x88x92W"sgr"(120))/(W"sgr"(60)+W"sgr"(120)xe2x88x922W"sgr"(0))}xc2x1xcfx80]xe2x80x83xe2x80x83(4)
Thus, at step 208, the expression of Eq. (4) is computed for each pixel based on the values of W"sgr"(0), W"sgr"(60), and W"sgr"(120) that were computed at step 206. Of the two solutions to equation (4), the direction xcex8max is then selected as the solution that yields the larger magnitude for W"sgr"(xcex8) at that pixel, denoted W"sgr"(xcex8max). Thus, at step 208, an array xcex8max(i,j) is formed that constitutes the direction image corresponding to the digital mammogram image I. As an outcome of this process, a corresponding two-dimensional array of line intensities corresponding to the maximum direction xcex8max at each pixel is formed, denoted as the line intensity function W"sgr"(xcex8max).
At step 210, a line image L(i,j) is formed using information derived from the line intensity function W"sgr"(xcex8max) that was inherently generated during step 208. The array L(i,j) is formed from W"sgr"(xcex8max) using known methods such as a simple thresholding process or a modified thresholding process based on a histogram of W"sgr"(xcex8max). With the completion of the line image array L(i,j) and the direction image array xcex8max(i,j), the line detection process is complete.
Optionally, in the Karssemeijer algorithm a plurality of spatial scale values "sgr"1, "sgr"2, . . . , "sgr"n may be selected at step 202. The steps 204-210 are then separately carried out for each of the spatial scale values "sgr"1, "sgr"2, . . . , "sgr"n. For a given pixel (i,j), the value of xcex8max(i,j) is selected to correspond to the largest value among W"sgr"1(xcex8max1), W"sgr"2(xcex8max2), . . . , W"sgr"n(xcex8maxn). The line image L(i,j) is formed by thresholding an array corresponding to largest value among W"sgr"1(xcex8max1), W"sgr"2(xcex8max2), . . . , W"sgr"n(xcex8maxn) at each pixel.
Although it is generally more computationally efficient than the filter-bank method, the prior art Karssemeijer algorithm has computational disadvantages. In particular, for a given spatial scale parameter "sgr", the Karssemeijer algorithm requires three separate convolutions of Nkxc3x97Nk kernels with the Mxc3x97N digital mammogram image I. Each convolution, in turn, requires approximately Mxc2x7Nxc2x7(Nk)2 multiplication and addition operations, which becomes computationally expensive as the kernel size Nk, which is proportional to the spatial scale parameter "sgr", grows. Thus, for a constant digital mammogram image size, the computational intensity of the Karssemeijer algorithm generally grows according to the square of the scale of interest.
Accordingly, it would be desirable to provide a line detection algorithm for use in a medical imaging system that is less computationally intensive, and therefore faster, than the above prior art algorithms.
It would further be desirable to provide a line detection algorithm for use in a medical imaging system that is capable of operating at multiple spatial scales for detecting lines of varying widths.
It would be even further desirable to provide a line detection algorithm for use in a medical imaging system in which, as the scale of interest grows, the computational intensity grows at a rate less than the rate of growth of the square of the scale of interest.
These and other objects are provided for by a method and apparatus for detecting lines in a medical imaging system by filtering the digital image with a single-peaked filter, convolving the resultant array with second order difference operators oriented along the horizontal, vertical, and diagonal axes, and computing direction image arrays and line image arrays as direct scalar functions of the results of the second order difference operations. Advantageously, it has been found that line detection based on the use of four line operator functions can actually require fewer computations than line detection based on the use of three line operator functions, if the four line operator functions correspond to the special orientations of 0, 45, 90, and 135 degrees. Stated another way, it has been found that the number of required computations is significantly reduced where the aspect ratio of the second order difference operators corresponds to the angular distribution of the line operator functions. Thus, where the second order difference operators are square kernels, having an aspect ratio of unity, the preferred directions of four line operator functions is at 0, 45, 90, and 135 degrees.
In a preferred embodiment, a spatial scale parameter is selected that corresponds to a desired range of line widths for detection. The digital image is then filtered with a single-peaked filter having a size related to the spatial scale parameter, to produce a filtered image array. The filtered image array is separately convolved with second order difference operators at 0, 45, 90, and 135 degrees. The direction image array and the line image array are then computed at each pixel as scalar functions of the elements of the arrays resulting from these convolutions. Because of the special symmetries involved, the second order difference operators may be 3xc3x973 kernels. Moreover, the number of computations associated with the second order difference operations may be achieved with simple register shifts, additions, and subtractions, yielding an overall line detection process that is significantly less computationally intensive than prior art algorithms.
In another preferred embodiment, the digital image is first convolved with a separable single-peaked filter kernel, such as a Gaussian. Because a separable function may be expressed as the convolution of a first one dimensional kernel and a second one dimensional kernel, the convolution with the separable single-peaked filter kernel is achieved by successive convolutions with a first one dimensional kernel and a second one dimensional kernel, which significantly reduces computation time in generating the filtered image array. The filtered image array is then convolved with three 3xc3x973 second order difference operators, the first such operator comprising the difference between a horizontal second order difference operator and a vertical difference operator, the second such operator comprising the difference between a first diagonal second order difference operator and a second diagonal second order difference operator, and the third such operator being a Laplacian operator. Because of the special symmetries associated with the selection of line operator functions at 0, 45, 90, and 135 degrees, the direction image array and the line image array are then computed at each pixel as even simpler scalar functions of the elements of the arrays resulting from the three convolutions.
Thus, line detection algorithms in accordance with the preferred embodiments are capable of generating line and direction images using significantly fewer computations than prior art algorithms by taking advantage of the separability of Gaussians and other symmetric filter kernels, while also taking advantage of discovered computational simplifications that result from the consideration of four line operator functions oriented in the horizontal, vertical, and diagonal directions.